Getting started

Preliminaries

Install and load ciftiTools:

if (!require("ciftiTools",  quietly=TRUE)){ 
  install.packages("ciftiTools")
} else if (packageVersion("ciftiTools") < "0.9.0") {
  update.packages("ciftiTools")
}

library(ciftiTools)

Download the Connectome Workbench here: https://www.humanconnectome.org/software/get-connectome-workbench

Tell ciftiTools where to find the Connectome Workbench:

# Replace '~/Desktop/workbench' with the actual path on your computer.
#   If successful, the path to the Workbench executable will be printed.
ciftiTools.setOption("wb_path", "~/Desktop/workbench")
## Using this Workbench path: '/Users/ddpham/Desktop/workbench/bin_macosx64/wb_command'.

Setup rgl plots for knitting.

library(rgl)
## Warning: package 'rgl' was built under R version 4.1.2
rgl::setupKnitr()

# Sometimes the first OpenGL window does not render properly.
rgl::rgl.open(); rgl::rgl.close()

First steps

Read in our first CIFTI file!

cii_fname <- file.path("data", "MyelinAndCorrThickness.dscalar.nii")
xii <- read_cifti(cii_fname)
xii # same as `summary(xii)`
## ====CIFTI METADATA===================
## Intent:           3006 (dscalar)
## - names           "MyelinMap_BC_decurv", "corrThickness"
## Measurements:     2 columns
## 
## ====BRAIN STRUCTURES=================
## - left cortex     30424 data vertices
##                   2068 medial wall vertices (32492 total)
## 
## - right cortex    30527 data vertices
##                   1965 medial wall vertices (32492 total)

Plot it:

# same as view_xifti(xii) and view_xifti_surface(xii)
plot(xii, zlim=c(1,2))
Myelin

Myelin

View a list of commonly-used functions:

if(interactive()) { help(ciftiTools) }

View details for a specific function:

if(interactive()) { help(view_xifti_surface) }

Task #1

a. Read & plot the fMRI data.

Use read_cifti to read in the CIFTI file named cii_fname. Set the brainstructures argument to include all brain structures.

cii_fname <- file.path("data/toy_fMRI.dtseries.nii")
# [Use `read_cifti` here]
xii <- read_cifti(cii_fname, brainstructures="all")
xii
## ====CIFTI METADATA===================
## Intent:           3002 (dtseries)
## - time step       1 (seconds)
## - time start      0
## Measurements:     500 columns
## 
## ====BRAIN STRUCTURES=================
## - left cortex     3670 data vertices
##                   332 medial wall vertices (4002 total)
## 
## - right cortex    3674 data vertices
##                   328 medial wall vertices (4002 total)
## 
## - subcortex       31870 data voxels
##                   subcortical structures and number of voxels in each:
##                     Cortex-L (0), Cortex-R (0),
##                     Accumbens-L (135), Accumbens-R (140),
##                     Amygdala-L (315), Amygdala-R (332),
##                     Brain Stem (3472),
##                     Caudate-L (728), Caudate-R (755),
##                     Cerebellum-L (8709), Cerebellum-R (9144),
##                     Diencephalon-L (706), Diencephalon-R (712),
##                     Hippocampus-L (764), Hippocampus-R (795),
##                     Pallidum-L (297), Pallidum-R (260),
##                     Putamen-L (1060), Putamen-R (1010),
##                     Thalamus-L (1288), Thalamus-R (1248).

We will add very inflated surfaces prior to plotting.

xii <- add_surf(
  xii, 
  load_surf("left", "very inflated"), 
  load_surf("right", "very inflated")
)

Use view_xifti_surface to plot the first timepoint.

Customize your plot by specifying some arguments! Here’s a few ideas:

  • Set title
  • Set zlim (min, max)
  • Change colors to the name of a sequential palette from RColorBrewer. To see a list, use RColorBrewer::display.brewer.all().
# [Use `view_xifti_surface` here]
view_xifti_surface(
  xii, title="T=1", 
  zlim=c(5000, 10000), colors="PuRd"
)
First timepoint

First timepoint

b. Compute & plot the fALFF

Here’s a function to calculate the fALFF of a single time series vector.

#' Adapted from `ANTsR::alffmap`
falff <- function(x, TR=.72, f_low=.01, f_high=.1){
  # Subtract linear trend
  x <- lm(x~seq(length(x)))$residuals
  
  # Compute Fourier decomposition.
  y <- stats::spec.pgram(
    stats::ts(x, frequency=1/TR), 
    taper=0, plot=FALSE
  )
  freq_mask <- (y$freq > f_low) & (y$freq < f_high)
  
  # Compute proportion within frequency range
  sum(y$spec[freq_mask]) / sum(y$spec)
}

Use apply_xifti to apply the falff function to each row of the data matrix in xii.

# [Use `apply_xifti` here.]
fxii <- apply_xifti(xii, 1, falff)

Plot the results using a color scale going from 0 (white) to 1 (blue), and save to a PNG file. For help on finding the relevant arguments, read help(view_xifti_surface).

# [Use `view_xifti_surface` here.]
view_xifti_surface(fxii, colors=c("white", "blue"), zlim=c(0, 1))
fALFF

fALFF 

# [Use `view_xifti_volume` here.]
view_xifti_volume(fxii, colors=c("white", "blue"), zlim=c(0, 1))
fALFF, subcortex

fALFF, subcortex 

view_xifti_surface(fxii, colors="Blues", zlim=c(0, 1))
fALFF, alternate colors

fALFF, alternate colors

Task #2

Read the fMRI data and parcellation.

cii_fname <- file.path("data/MSC03_rest07.4k.dtseries.nii")
xii <- read_cifti(cii_fname)
dim(xii)
## [1] 7344  818
parc <- read_cifti("data/Kong2022_100Parc.dlabel.nii", resamp_res=4000)
dim(parc)
## [1] 8004    1
# rownames(parc$meta$cifti$labels[[1]])
table(c(as.matrix(parc))==0)
## 
## FALSE  TRUE 
##  7219   785

Incompatible ROIs: we need to replace the medial wall mask with NA values.

xii <- move_from_mwall(xii)

Pick our seed and plot it.

# View(parc$meta$cifti$labels$parcels["Key"])
plot(parc, title="Parcellation")
Parcellation

Parcellation 

the_seed <- 3 # sample(seq(100), 1) 
parc_masked <- transform_xifti(parc, function(q){ifelse(q==the_seed, q, 0)})
plot(parc_masked, borders="black", title=paste0("Parcel #", the_seed))
Selected parcel

Selected parcel

Now that we’ve identified the seed, let’s do this!

a. Scrub the fMRI data.

Use select_xifti to mask out high-motion timepoints. (Only keep the columns for which the corresponding entry of xii_tmask is TRUE.)

xii_tmask <- read.csv("data/MSC03_rest07_tmask.txt")
xii_tmask <- which(as.logical(xii_tmask[,1]))
# [Use `select_xifti` here.] 
xii <- select_xifti(xii, xii_tmask)
dim(xii)
## [1] 8004  475

b. Compute the seed correlations.

The below code chunk is already filled out. Here’s what it does:

  • seed_mask: Convert the parcellation "xifti" to a numeric vector of parcel indices for each vertex, and then convert that to a binary vector indicating whether the parcel is the seed parcel.
  • parc_mean_ts: Convert the resting-state fMRI "xifti" to a matrix, subset the rows to select only vertices inside the parcel, and then take the mean across space.
  • xii_sc: Compute correlation between each vertex time series, and the parcel mean time series.
seed_mask <- c(as.matrix(parc)) == the_seed
parc_mean_ts <- colMeans(as.matrix(xii)[seed_mask,,drop=FALSE])
xii_sc <- cor(t(as.matrix(xii)), parc_mean_ts)
xii_sc <- newdata_xifti(xii, xii_sc)

Now we will convert the result back to a "xifti" with newdata_xifti, and then format is as a “dscalar” since the data no longer represents a time series.

xii_sc <- cor(t(as.matrix(xii)), parc_mean_ts)
# [Use `newdata_xifti` here.]
xii_sc <- newdata_xifti(xii, xii_sc)
# [Use `convert_xifti` here.]
xii_sc <- convert_xifti(xii_sc, "dscalar")

Write out the result with write_xifti.

seed_fname <- paste0("seed_cor_", the_seed, ".dscalar.nii")
# [Use `write_xifti` here.]
write_xifti(xii_sc, seed_fname)
## Writing left cortex.
## Writing right cortex.
## Creating CIFTI file from separated components.

Finally, plot the result!

view_xifti_surface(xii_sc, zlim=c(-.7, .7))
Seed correlation

Seed correlation

Task #3

Install the latest version of BayesfMRI from GitHub.

if (!requireNamespace("BayesfMRI", quietly=TRUE) || packageVersion("BayesfMRI") < "0.2.0") {
  if (!requireNamespace("devtools", quietly=TRUE)) { install.packages("devtools") }
  devtools::install_github("mandymejia/BayesfMRI", "2.0")
}
library(BayesfMRI)

Read in and format the event onsets file. Use the “left foot” task.

onsets <- paste0("data/MSC03_motor07_run0", seq(2), "_events.tsv")
onsets <- lapply(as.list(onsets), read.csv, sep="\t")
onsets <- lapply(onsets, function(q){list(as.matrix(q[q[,3]=="LFoot",seq(2)]))})

Run classical GLM for the right and left cortex separately.

bglmR <- BayesGLM_cifti(
  paste0("data/MSC03_motor07_run0", seq(2), ".4k.dtseries.nii"),
  brainstructures="right", surfR = load_surf("right", resamp_res=4000),
  onsets=onsets, ar_order=0, EM=FALSE, Bayes=FALSE, resamp_res=NULL, TR=2.2
)
## 
##  SETTING UP DATA 
##  .. reading in data for session 1
##  .. reading in data for session 2
##  MAKING DESIGN MATRICES 
##  RUNNING MODELS 
## 
##  .. RIGHT CORTEX ANALYSIS 
##   328 locations removed due to NA or NaN values in at least one scan.
##  PUTTING RESULTS IN CIFTI FORMAT 
## 
##  DONE!
bglmL <- BayesGLM_cifti(
  paste0("data/MSC03_motor07_run0", seq(2), ".4k.dtseries.nii"),
  brainstructures="left", surfL = load_surf("left", resamp_res=4000),
  onsets=onsets, ar_order=0, EM=FALSE, Bayes=FALSE, resamp_res=NULL, TR=2.2
)
## 
##  SETTING UP DATA 
##  .. reading in data for session 1
##  .. reading in data for session 2
##  MAKING DESIGN MATRICES 
##  RUNNING MODELS 
## 
##  .. LEFT CORTEX ANALYSIS 
##   332 locations removed due to NA or NaN values in at least one scan.
##  PUTTING RESULTS IN CIFTI FORMAT 
## 
##  DONE!

a. Write combine_and_plot_GLM

This function should combine classical, multi-session GLM results from the left and right hemispheres, and then plot the result. To make it work, add a line that uses combine_xifti to merge the left and right hemisphere "xifti" objects.

combine_and_plot_GLM <- function(
  left_GLM, right_GLM, 
  what=c("estimates", "SE_estimates", "resids"),
  ...){
  
  # Argument checks.
  stopifnot(is.xifti(left_GLM$betas_classical$session_avg))
  stopifnot(is.xifti(right_GLM$betas_classical$session_avg))
  what <- match.arg(what, c("estimates", "SE_estimates", "resids"))
  
  # Get separate `"xifti"` objects.
  if (what=="estimates"){
    left_GLM <- left_GLM$betas_classical$session_avg
    right_GLM <- right_GLM$betas_classical$session_avg
    
  } else if (what=="SE_estimates") {
    # Left
    left_mask <- left_GLM$GLMs_classical$cortexL$session_avg$mask
    x <- left_GLM$GLMs_classical$cortexL$session_avg$SE_estimates
    left_GLM <- newdata_xifti(
      left_GLM$betas_classical$session_avg,
      x[left_mask,,drop=FALSE]
    )
    # Right
    right_mask <- right_GLM$GLMs_classical$cortexR$session_avg$mask
    x <- right_GLM$GLMs_classical$cortexR$session_avg$SE_estimates
    right_GLM <- newdata_xifti(
      right_GLM$betas_classical$session_avg,
      x[right_mask,,drop=FALSE]
    )
    
  } else if (what=="resids") {
    left_GLM <- newdata_xifti(
      left_GLM$betas_classical$session_avg,
      left_GLM$GLMs_classical$cortexL$session_avg$resids
    )
    right_GLM <- newdata_xifti(
      right_GLM$betas_classical$session_avg,
      right_GLM$GLMs_classical$cortexR$session_avg$resids
    )
  } else { stop() }
  
  # [Use `combine_xifti` here.]
  both_GLM <- combine_xifti(left_GLM, right_GLM)
  
  plot(both_GLM, ...)
}

Try it out!

combine_and_plot_GLM(bglmL, bglmR, "SE")
## `zlim` not provided: using color range 0 - 0.598 (data limits: 0.0725 - 3.58).

Lastly, let’s check out an example of using view_comp. We will plot the GLM result alongside its activations calculated by statistical inference.

# Write and save the GLM result plot.
p1 <- plot(bglmR, fname="glmR_vals", title="GLM result", together="leg")
# Write and save the activations plot.
p2 <- plot(id_activations_cifti(bglmR)$activations_xifti, fname="glmR_act", title="Thresholded activation", together="leg")
# Composite
view_comp(c(p1, p2), nrow=1)
Left foot task activation (right cortex)

Left foot task activation (right cortex)

Appendix

Bonus Task #1

Read in cortex CIFTI data without using read_xifti! To work around this challenge, use separate_cifti to split the CIFTI file into GIFTI files, and then read those in with read_xifti2. Check your work by comparing the result to read_xifti. Extra bonus if you can find a way to obtain the subcortical data too.

xii1 <- read_xifti("data/MSC03_motor07_run01.4k.dtseries.nii")

xii2_sep <- separate_cifti("data/MSC03_motor07_run01.4k.dtseries.nii")
xii2 <- read_xifti2(
  xii2_sep["cortexL"], xii2_sep["ROIcortexL"],
  xii2_sep["cortexR"], xii2_sep["ROIcortexR"]
)

stopifnot(max(as.matrix(xii1-xii2)) == 0)
# Extra bonus
xii1 <- read_xifti(
  "data/MSC03_motor07_run01.4k.dtseries.nii", 
  brainstructures="all"
)

xii2_sep <- separate_cifti(
  "data/MSC03_motor07_run01.4k.dtseries.nii",
  brainstructures="all"
)
xii2_ctx <- read_xifti2(
  xii2_sep["cortexL"], xii2_sep["ROIcortexL"],
  xii2_sep["cortexR"], xii2_sep["ROIcortexR"]
)
xii2_sub <- as.xifti(
  subcortVol = RNifti::readNifti(xii2_sep["subcortVol"]),
  subcortLabs = RNifti::readNifti(xii2_sep["subcortLabs"]),
  subcortMask = RNifti::readNifti(xii2_sep["ROIsubcortVol"])
)
xii2 <- combine_xifti(xii2_ctx, xii2_sub)
stopifnot(max(as.matrix(xii1-xii2)) == 0)

Making plots for publications

ciftiTools can create high-quality plot suitable for publication. Here are a few tips:

  • Save PNG files using fname rather than taking screenshots of the OpenGL window.
  • Use width and height to increase the plot size. The limits are equal to your screen’s resolution.
  • view_comp can create composite images, but it’s also possible to put together plots by importing them into an image-editing software like GIMP, or to create subplots within LaTeX. For this same reason, you may want to omit title and instead add text above or below the plot later. You can also use legend_embed=FALSE to composite the legend yourself.
  • When plotting dlabel files where colors are repeated for different keys, you can use borders=TRUE to delineate them better on cortex plots.
  • In Photoshop, the “Object Selection” tool works great with the brain meshes! You can also try using the material=list(smooth=FALSE, lit=FALSE) to avoid smoothing colors and shading, allowing precise selection of region by color/value for further editing.

Please see the Appendix C in the ciftiTools paper for guidance on embedding ciftiTools plots in R Markdown documents.